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Abstract 

. In this paper we construct a third order method for solving additively split autonomous stiff 



systems of ordinary differential equations. The constructed additive method is L-stable with respect 
to the implicit part and allows to use an arbitrary approximation of the Jacobian matrix. In opposite 
' to our previous paper [8], the fourth stage is explicit. So, the constructed method also has a good 

, stability properties because of L-stability of the intermediate numerical formulas in the fourth 

stage, but has a lower computational costs per step. Automatic stepsize selection based on local 
error and stability control are performed. The estimations for error and stability control have been 
obtained without significant additional computational costs. Numerical experiments show reliability 
" <^ , and efficiency of the implemented integration algorithm. 

m 

■ 1 Introduction 

cn 
(N 

' difference or finite element methods results in the Cauchy problem for the system of ordinary differential 

O ■ equations with an additively split right hand side function of the form: 



Spatial discretization of continuum mechanics problems in partial differential equations by finite 



y' = ^{t,y) + 9{t,y), y(io)=yo, to<t<tk, 

where (p{t, y) is a non-symmetrical term obtained from discretization of the first-order differential 
operator, g{t, y) is a symmetrical term obtained from discretization of the second-order differential 
operator, t is a independent variable. It is assumed that in the problem the vector-function (7 is a stiff 
term and is a non-stiff term. 

Explicit Runge-Kutta methods have a bounded stability region and are suitable for non-stiff and 
mildly stiff problems only. L-stable methods are usually used for solving stiff problems. In the case of 
large-scale problems overall computational costs of L-stable methods are almost completely dominated 
by evaluations and inversions of the Jacobian matrix of a right hand side vector function. Overall 
computational costs can be significantly reduced by re-using the same Jacobian matrix over several 
integration steps (freezing the Jacobian). 

Freezing the Jacobian in iterative methods has effect on convergence speed of an iterative process 
only and doesn't lead to loss of accuracy. So, this approach is extensively used for implementation 
of these methods. For Rosenbrock type methods and their modifications [4] an approximation of the 
Jacobian matrix can lead to decreasing a consistency order. 

The system y' = f{t,y) can be written in the form y' = [f{t,y) — By] + By, where B is 
some approximation of the Jacobian matrix. Assume that stiffness is fully concentrated in the term 
g{t,y) = By, then the expression (p{t,y) = f{t,y) — By can be interpreted as the non-stiff term [2,7]. 
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If the Cauchy problem is considered in the form y' = [f{t, y) — By] + By under construction of additive 
methods, then an arbitrary approximation of the Jacobian matrix can be used without decreasing the 
order of these methods. Additive methods constructed in this way allow both analytical and numerical 
computations of the Jacobian matrix. Note that the approximation of the Jacobian by a diagonal 
matrix is suitable for some mildly stiff problems. 

In this paper we construct a six-stage third order additive method that allows to use different 
kinds of approximation of the Jacobian matrix. In opposite to our previous paper [8], the fourth stage 
is explicit. The constructed method also has a good stability properties because of L-stability of the 
intermediate numerical formulas (with respect to the implicit part) in the fourth stage, but has a lower 
computational costs per step. The estimation of the error has been obtained on the base of an embedded 
additive method without any additional computational costs. The estimation of the maximum absolute 
eigenvalue of the Jacobian matrix has been obtained by a power method using only two additional 
computations of </?(y). Hence, additional computational costs will be negligible, especially for large- 
scale problems. These estimations are used for error and stability control correspondingly. Numerical 
experiments are performed showing the reliability and efficiency of the constructed method. 



2 A numerical scheme for autonomous problems 

Consider the Cauchy problem for an autonomous system of ordinary differential equations 

y' = (p{y) + g{y), y[to) = yo, to < t < tj,, (1) 

where y, 99 and g are A'^- dimensional smooth vector-functions, t is an independent variable. In the 
following, we assume that g is a stiff term and </? is a non-stiff term. Consider a six-stage numerical 
scheme for solving (1): 

6 

Vn+l =yn + '^Piki , 
i=l 

ki = hip{yn), 
Dnk2 = h[ip{yn) +g{yn)], 

Dnkz = k2, 

3 3 

ki = h(p{yn + ^ Pijkj) + hg{yn + ^ a^jkj), (2) 
Dnk^ = ki + 7/03, 

5 

ka = hifiyn + ^^6i%), 
i=i 

where -D„ = E — ahg'^, g'^ = dg{yn)/dy is the Jacobian matrix of the function g{y), E is the identity 
matrix, fcj, 1 < i < 6, are stages, a,pi,a4j, l34^j, Pqj,^ are coefficients that have effect on accuracy and 
stability properties of the scheme (2). 



3 The third order conditions 

The Taylor series expansion of the approximate solution up to terms in has the form 



Vn+l =yn+ {Pl + P2 + P3 + P4 + (7 + 1)P5 + Po) hip + {p2 + Ps + P4 + (t + 1)P5) hg+ 
+ ((/341 + ^42 + Paz){pa + P5) + (/?61 + /362 + ^63 + PdA + (t + 1)/?65)P6) /lVV+ 
+ ((/?42 + P43){P4 + P5) + {Pe2 + Pes + 064 + {l + l)P65)Pe)h^V^' 9+ 
+ [a{p2 + 2p3 + (37 + 1)^5) + (0:41 + a42 + Q;43)(P4 +P5)]h^9'v>+ 
+ [a{p2 + 2p3 + (37 + 1)^5) + (a42 + q;43)(P4 +P5)]h^g'g+ 

+ 0.5[(/341+/342+/?43)^(P4+P5) + (/?61+/962+/?63 + /364 + (7+l)/?65)^P6] /iVV^ + 
+ 0.5[(/?42 + (343)\P4 + P5) + {P62 + p63 + ^64 + (7 + l)/?65)'p6] + 
+ [{P42 + P43){P41 + /342 + ^43)^4 + P5) + {PgI + ^62 + /363 + ^64+ 
+ (7 + l)/965) (/?62 + P63 + /?64 + (7 + l)/?65)p6] h^^"vg+ 

+ (/?41 + /342 + /343)(/964 + f3e5)P6h^^'% + (/342 + /343)(/364 + Mp^h"^ ^'"^ 9+ 
+ [a((/?42 + 2/343) {P4+Pr>) + {I3%2 + 2/363 + (37 + 1)/365)P6) + 

+ (q;4i + a42 + 043) (/364 + /365)P6 h^^'g'^ + a((/342 + 2/343) (P4 + Pb) + 

+ (^62 + 2/363 + (37 + 1)^65)P6) + ("42 + "43) {P&4 + I3m)p&\ h^v'g'g+ 

+ 0.5(a4i + 0:42 + 0:43)^ (P4 + P5)h^g"<P^ + 0.5(a42 + ^43)^(^4 + Pb)h^g"g'^+ 

+ (a41 + 042 + Q!43)(Q;42 + 043) (P4 + Pb)h^g"'Pg + o(/341 + /342 + /343)P5^^5W+ 
+ a(/342 + P43)P5h^g''p'g + a[a{p2 + Sps + (67 + 1)^5) + (042 + 2q;43)p4+ 
+ (Q!41 + 2q;42 + 30:43)^5] h^g''^'P + a[a{p2 + ^P3 + (67 + 1)P5) + ("42 + 20:43)^4+ 
+ (2^42 + 3a43)P5] h^g'^g + 



where the corresponding elementary differentials are evaluated at i/n- 

The Taylor series expansion of the exact solution up to third order terms is 

/j2 ^3 
y(Wi) = y(tn) + h{ip + g) + —{ip'ip + ip'g + g'ip + g'g) + —iv"<p^+ 

+ ip"g'^ + 2ip"^g + + ¥''^5 + v'g'v + g' g + ^'V^ + (3) 

+ 2g"^g + 5VV + g'^'g + + g'''g) + o(^^), 

where the corresponding elementary differentials are evaluated at y{tn). 

Comparing the successive terms in the Taylor series expansion of the approximate and the exact 
solutions up to third order terms under the assumption ?/„ = y{tn) we have the system of nonlinear 
algebraic equations. Its solving results in the relation /34i = 041 = /36i = and the third order 



conditions of the scheme (2) take form: 

P2 + P3 + P4 + (7 + 1)P5 = 1, 
a(^42 + f3i3)p5 = 1/6, 

(/342+/343)(/364+/?65)P6 = 1/6, 

(/342 + /343)(P4 + P5) + [/?62 + /363 + /364 + (t + 1)/365]P6 = 0.5, (4) 
(/342 + /343)^(P4 + Ps) + [P62 + p63 + PgA + (t + l)/365] V = 1/3, 

a(^42 + 2/343) (P4 + Ps) + [a (/362 + 2/?63 + (37 + 1)/365) + ("42 + "43) (/?64 +/365)] P6 = 1/6, 
(042 + a43)^(P4 +P5) = 1/3, 

a[p2 + 2p3 + (37 + l)p5) + (a42 + a43)(P4 +P5) = 0.5, 

a[a{p2 + 3p3 + (67 + 1)^5) + (042 + 2a43)p4 + (2042 + 3043)^5] = 1/6, 

an = /?4i = Pqi =0, pi = -pQ. 



4 Stability analysis 

The hnear stabihty analysis of the additive scheme (2) is based on the scalar model equation 

y' = \iy + hy, y(0) =yo, t> 0, Re(Ai) < 0, Re(A2) < 0, | Re(Ai)j « | Re(A2)|, (5) 

where the free parameters Ai, A2 can be interpreted as some eigenvalues of the Jacobian matrices of 
the functions f (the non-stiff term) and g (the stiff term) correspondingly. 

Application of the scheme (2) for numerical solving the equation (5) yields 

Vn+i = R{x,z)yn, 

where x = Ai/i, z = A2/1 and R{x,z) is a stability function (its analytical expression is omitted here 
for brevity). 

The necessary condition of L-stability of the additive scheme (2) with respect to the stiff term 
has the form: 

lim R(x, z) = 0. 
2— >— 00 

It is satisfied if the following conditions hold: 

0(42 = a, (3^2 = 0, 

o?{pi + pq) - a(/?62 + /?64)P6 + a^sPeiPe = 0, (6) 
- a{p2 + Pa) + 043^4 = 0. 

Solving the system (4), (6). In the following, we assume that ^^j—i ot^j — 1 , (3q2 — o. The 

3 

first relation ensures that g{yn + Yl ^'ij^^j) approximate g{y{tn+i)) in the fourth stage and the other 

one improve stability properties of the intermediate numerical formula. 
Let us denote 

(5i = Pei + (365 , P2 = P62 + P63 + P6A + (7 + l)/?65 , 

P3 = a{P(i2 + 2/363 + (37 + l)/?65) + PdA + P65, U = 0(7 - 1) + 1. 



Then after obvious simplifications the system (4), (6) takes the form 



P2 + P3 + P4 + (7 + 1)P5 = 1, 

aPisPb = 1/6, 

PAsiPsi + Pe5)P6 = 1/6, 

/?43(P4 + P5) + {P62 + (im + /364 + (t + 1)/365)P6 = 0.5, 
PlziPi + P5) + {062 + (363 + /?64 + (t + 1)/365)^P6 = 1/3, 

2a/343(P4 + Ps) + (a(/?62 + 2/363 + (37 + 1)P65) + /?64 + P65)P6 = 1/6, (7) 
P4+P5 = 1/3, 

ap2 + 2a^»3 +P4 + (a(37 + 1) + l)p5 = 0.5, 
a{ap2 + 3ap3 + (2 - a)p4 + 3{2a-y + 1)^5) = 1/6, 

— ap2 + (1 — 2a)p4 = 0, 
a4i = /?4i = /342 = Pel = 0, 042 = 062 = a, 043 = 1 - a, /364 = a^/ (1 - 2a), pi = -pQ. 

Multiplying the first equation of (7) by 2 and subtracting the result from the eight one we obtain 
—ap2 + (1 — 2a)p4 + np5 = 0.5 — 2a. It follows from here and the tenth equation of (7) that 

P5 = 0.5(20^ -4a + l)/u. (8) 

We shall try to obtain an equation for a. For this purpose we divide the ninth equation by a and 
subtract the eighth one from the result. As the result we obtain: aps = (a — l)p4 + (a — 807 — 2)p5 + 
(1 — 3a) /{6a). Substituting this relation to the eighth equation we have —ap2 = {2a — l)p4 — 3up5 + 
(2 — 9a) /{6a). It follows from here and the tenth equation of (7) that (6a^ — 9a + 2)/(6a) — 3np5 = 0. 
Substituting this relation to (8) we obtain the following equation for a: 

6a^ - ISa^ + 9a - 1 = 0. (9) 
Then from the second equation of (7) we have 

P43 = l/(6ap5), (10) 
It follows from (10) and the third equation of (7) that 

Pi = ap5/p6- (11) 
Prom (11) and the notation Pi = Pq4 + /Jes we obtain 

Peb = ap5/pe - Pm, (12) 
from (10) and the fourth and the seventh equations of (7), we have 

p2 = l/{2p6)-l/{18ap5P6)- (13) 

It follows from sixth and seventh equations of (7) and (10) that P3 = 1/{6pq) — 1/{9p5Pq). Prom (12), (13) 
and the relations P2 = P62 + Pes + Pei + (7 + l)/365, Pe2 = a, we obtain 

Pes = 0.5(1 - 2a(7 + 1)P5)/P6 - l/(18ap5P6) - a + 7/364- (14) 



It follows from Ps = a(/362 + 2/363 + (37 + l)/365) + Pga + /^es, (12), (14), the eleventh equations of (7) 

that: ^(0^(7 — l) + o)p5 + a — 1/6^ /pq = (0^(7 — l) + a^) /(I — 2a). Using the notations introduced above 

the last equation can be written in the form: (aups + a — 1/6) /pg = a^ti/(l — 2a). Substituting (8) 
to this equation, we have (6a^ — 12a^ + 9a — 1)/{6pq) = a^u/(l — 2a). Next, using (9) we obtain the 
following expression for pQ : 

P6 = (1 - 2a)/n. (15) 

It follows from the obtained results, the seventh equation of (7) and (8) that p4 = (— 6a^ + 2a(7 + 
5) — 1)/ (6ii). Substituting the last relation to the tenth equation of (7) we obtain the expression for p2 : 
IIoflCTaBjiHH nocjiefluee paseucTBO b ^ecHToe ypaBHeHue (7), BbipasuM p2 '■ P2 = (6a^(7 + 1) — 4a^(7 + 
5) + 2a(7 + 6) — 1) /{Qau). It follows from here and from (9) that 

P2 = (2a2(77-l)-a(77-3)+7)/(6a'u). (16) 

Substituting the seventh equation of (7), (8) and (16) to the first one we express ps : P3 = (— 6a^7 + 
2a2(7-l) + a(47 + l)-7)/(6aw). Now, using (9), we have ps = (-20^(87 + 1) + 0(187 + 1)- 27) /(6au). 
Prom (8), (11) and (15), we have jSi = 0.5a(2a^ — 4a + 1)/(1 — 2a). It follows from here and from (9) 
that /?! = (6a^ — 6a + l)/(6 — 12a). Prom the last relation, the eleventh equation of (7) Pq4 = a^/(l — 2a) 
and /3i = (3q4 + /Jgs we obtain /Jgs = —(6a — l)/(6 — 12a). From (8), (9) and (10) we have 

/543 = n/(6a2-6a + l). (17) 

Substituting (3g4, = aV(l-2a), (8) and (15) to (14) and using (9), we obtain (3^3 = {-2a^-/^ + 2{53a^ - 
35a + 4)7 + 142a2 - 86a + 9)/(6(-18a2 + 10a - 1)) . 

We shall try to obtain an equation for 7. From the fourth and the fifth equations of (7) we 
obtain (32 = 2(1 — (3l^)/{'i — 2/^43). It follows from here and from the fifth equations of (7) that 
P6 = (3 — 2/^43)^/(12 — 12Pl^). Substituting (17) to the last relation and using (9), we have 

_ 994a^ - 72a^7 + 4a^7(7 + 16) - 4a(7 + 143) + 67 

12 (I02a2 - a2(7 - 1)2 _ 2a(7 + 29) + 6) ' ^ ' 

Comparing (15) with (18), and using (9), we obtain the following equation for 7: 

4a^7^ - 12a(15a2 - 10a + 1)7^ + 3a(374a2 - 228a + 33)7 + 4(8130^ - 486a + 175/3) = 0. (19) 

Now, the coefficients of the L -stable third order scheme (2) can be computed by the following formulas: 



"41 


— (3ai — f3i2 — Pel — 0, 




"42 = /362 = a, 043 = 1 - a, 


Pi 


= -(l-2a)/^.. 




P2 = (2a2(77-l)-a(77-3)+7)/(6an), 


P3 


= (-2a2(87 + l) +a(137 + l) - 


27)/(6an), 


P4 = (-6a2 + 2a(7 + 5) - l)/(6u), 


P5 


= 0.5(2a2 -4a + l)/u, 




P6 = (1 - 2a)/n, 


Pi3 


= l/(6ap5). 




/3i = ap5/p6, (20) 


P2 


= 2(1 -/323)/(3- 2/343), 




Pu = aV(l - 2a), 


P65 


= /3l — /364, 




Pm = P2-a- P&A - (7 + l)/365- 


where u 


= a(7 — 1) + 1, the coefficients 


a and 7 is 


determined from the equations (9) and (19) 



correspondingly. 



The equation (9) has the following tree real roots: 

01 = 0.15898389998867, 02 = 0.43586652150845, 03 = 2.40514957850286. 

The numerical experiments show that the root 02 is the most suitable. The equation (19), in turn, has 
the following tree real roots under the condition = 02: 

^2, 1 = -4.51745281449726, 72, 2 = -2.49646456973997, 72, 3 = -1.02332630944762 . 

The numerical experiments show that the root 72,1. is the most suitable. Therefore computational 
results will be given for o = 02 and 7 = 72, i- 

The corresponding coefficients of the L -stable third order scheme (2) take the form 



o = +0.43586652150846, 7 

= +0.09130146290929, p2 

P4 = +0.20395977226114, ps 

a4i = 0, 042 

Pii = 0, (3i2 
Pel = 0, 

/3g4 = +1.48112677684356, /Jes 



-4.51745281449727, 

+0.49588787677190, ps 

+0.12937356107220, pe 

+0.43586652150846, 043 

0, Pa3 

+0.43586652150846, Pes 
-2.09874671679705 . 



+0.75521774748189, 
-0.09130146290929, 
+0.56413347849154, 
+2.95562753995095, 
-3.98487214709651, 



(21) 



5 Local error estimation 

For the error estimation we construct the embedded method of second order of the form: 

3 6 _ 

Vn+l, 2 = yn + '^ riki + ^ Tiki , 
i=l i=4 

ki = hip{yn), 

Dnk2 = h{ip{yn) + g{yn)) , 

Dnks = k2, (22) 

_ 3 

ki = hip{yn + ^ Pijkj), 

Dnk5 = /C4 + 7^3, 

_ 3 _ _ 

fee = hip{yn + ^ Pejkj + Pmh + P65h)- 
i=i 

where the coefficients r,, 1 < i < 6, should be determined, and parameters a, fS^j, f3Qj,'y are given 
by (20) or (21). Note that there is not hg{yn + Z]j=i ^ij^j) ™ the fourth stage as opposed to (2). 

The Taylor series expansion of the approximate solution computed by the scheme (22) up to terms 
in /i^ has the form 

Vn+i, 2 = yn+ [ri + r2 + r^+ri + {-f + l)r^ + tq) hip + (r2 + r-s + 'yr5)hg+ 
+ a(r2 + 2r3 + (37 + l)r5)h^g'ip + a(r2 + 2r3 + 3jrr,)h^g'g+ 

+ ((/341 + /?42 + P43){r4 + rs) + {Pel + Pm + /Ses + PeA + (7 + ^)Peb)re) 
+ {{Pi2 + Piz){ri + rs) + {Pe2 + Pm + lPe^)re)h^ip' g + 0{h^). 



where the elementary differentials are evaluated at yn- Comparing successive terms in the Taylor series 
expansion of the approximate and the exact solutions up to second order terms under the assumption 
Un = y{tn) we obtain the second order conditions of the scheme (22): 

1. ri + r2 + ra + r4 + (7 + l)r5 + rg = 1, 

2. r2 + rs + = 1, 

3. a(r2 + 2r3 + (37 + l)r5) =0.5, 

4. a{r2 + 2r3 + 37r5) = 0.5, (23) 

5. (/?4i + P42 + P43){r4 + r-s) + (/?6i + /?62 + Pas + /?64 + (7 + I)l365)r6 = 0.5, 

6. (/?42 + P43){r4 + rs) + (/?62 + ^63 + 1^65^6 = 0.5. 

Subtracting the third equation of (23) from the fourth one and using a / we obtain = 0. 
Subtracting the sixth equation from the fifth one and using the eleventh one of (7) and (3q4 + Pq^ ^ 
we have rg = 0. It follows from the second and the fourth equations of (23) that r2 = 0.5(4a — 
l)/a, ra = 0.5(1 — 2a)/a. Now, from the sixth and the first equation of (23) we have r4 = O.5/P43, 
n = -O.5//343. 

As the result we have all the coefficients of the L stable embedded method (22) of second order. 
For the coefficients (21) we obtain 

n = -0.16916881211910, r2 = +0.85285981986048 rg = +0.14714018013952, 

r4 = +0.16916881211910, = 0. rg = 0. 



The embedded method (22) doesn't require any additional computations of right hand side, 
evaluations and inversions of the Jacobian matrix, because of rg = and = 0. 
Let us denote the error estimation by 

IVn ~ ^n, 2I 

err„ = max — — ; — — — r-, 

i<i<N Atoli + Rtoli\yl\ 

where Atok and Rtok are the desired tolerances prescribed by the user. If err„ < 1, then the computed 
step is accepted, else the step is rejected and computations are repeated. When Rtok = 0, the absolute 
error is controlled on the i-th component of the solution with the desired tolerance Atoli. If Atoli = 
then the relative error is controlled on the i-th component with the tolerance Rtoli. 



6 Stability control and stepsize selection 



In the additive method (2) for solving (1) the non-stiff term (p is treated by the tree-stage exphcit Runge 
Kutta method (the exphcit part), and the stiff term g is treated by the L-stable (4,2)-method [9-11] 
(the imphcit part). In the general case there is no guarantee that the function 99(2/) = f{y) — By is the 
non-stiff term in reducing y' = f{y) to y' = [f{y) — By] + By. If some stiffness is in ip{y) = f{y) — By 
(i.e. stiffness leakage phenomenon occurs) then the additional stability control of the explicit part of 
the scheme (2) can increase efficiency of computations for many problems. In some cases it has no a 
significant effect on the efficiency of the integration algorithm because of the good stability properties 
of the scheme (2). Therefore the choice of using or not using the additional stability control of the 
explicit part is given to the end-user. 

We perform the stability control of the explicit part of the scheme (2) by analogy with [8]. For 
additive methods in opposite to explicit Runge Kutta methods it isn't possible to use previously 
computed stages because of peculiarity of the problem (1). Therefore instead of using the stages 
h, 1 < i < 6, of (2) we consider the additional stages di, ^2 of the form: 

di = hipiyn + a2iki), c?2 = hip{yn + asih + a32di). 

Denote ^p{y) = Ay + b, where A and b are matrix and vector with constant coefficients correspondingly, 
then we have 

ki = h{Ayn + b), di = ki + a2ihAki , d2 = ki + {aai + a32)hAki + a2ia32h^A^ki. 

Assuming 021 = 031 + 032 we obtain 

d2 — di = a2ia32h^A'^ki, di — ki = a2ihAki. 

The maximum absolute eigenvalue Vn = h\Xn max\ of the matrix hA can be approximated using the 
power method by the following formula: 

-1, 14-41 



1032 I max 



i<i<N \d\-k\\ 

then the stability control can be made by f„ < 2, where number 2 is an approximate length of the 
stability interval of the tree-stage explicit Runge Kutta method. 

In the general case this estimation is quite crude because of small number of iterations of the 
power method and the nonlinearity of the function ip{y). Therefore the stability control is used for 
limiting the stepsize growing only. 

Let the approximate solution y„ is computed with the stepsize hn- For the stepsize selection we 
use errn = 0{h^). The stepsize hacc predicted by accuracy we compute by the formula: hacc = Qihn, 
where qi is a root of the equation qferVn = 1- In view of Vn = 0{hn), the stepsize hst predicted by 
stability is computed by h^t = 92^n, where 52 is a root of the equation q2Vn = 2. Then the stepsize 
predicted by accuracy and stability is selected by the formula: 

hn+i = max[hn,inm{hacc,hst)]- 



The stability control of the explicit part of the scheme (2) requires, at each integration step, two 
additional computations of (p{y). These computational costs are negligible for large-scale problems, but 
if you are sure that all stiffness is in g{y) then you can take off stability control to save computational 
costs. 



7 Numerical experiments 



Further, the numerical code based on the additive method (2) (with error and stabihty control as 
well as with diagonal Jacobian approximation) is called AS0DE3 (the Additive Solver of Ordinary 
Differential Equations) . 

The test problems given below have been reduced to the form y' = {f{y) — By) + By. All numerical 
computations have been performed in double precision arithmetic on IBM PC Athlon(tm) XP 2000+ 
with the desired tolerances of the error Atol = Rtol = Tol = 10~™, m = 2,4. The scheme (2) is of 
third order, therefore it is unreasonable to do numerical computations with higher tolerance. 

The following four test examples are considered: 

Example 1 [17]. 

y[ = -0.013yi - lOOOyiya, 

y'^ = -2500^22/3, (24) 
y'^ = -0.013yi - lOOOyiya - 2500^2^3, 

t€[0,50], yi(0) = l, 2/2(0) = !, ^3(0) = 0, /lo = 2.9 • 10"^ 

Example 2 [13]. 

y[ = 77.27(^2 - yiy2 + yi- 8.375 • IQ-^y^), 

y'2 = {-y2 - 2/12/2 + 2/3)/77.27, (25) 
y^ = 0.161(yi-y3), 

t G [0, 300], yi(0) = 4, 1/2(0) = 1.1, ^3(0) = 4, /iq = 2 • 10-^. 



Example 3. 

y[ = -0.04yi + 0.0l2/22/3, 

y'2 = 400yi - 1002/2?/3 - 3OOO2/I, 

2/3 = 30yi, 

t G [0,40], yi(0) = l, y2(0) =y3(0) = 0, /iq = 10"^ 

Example 4. 

y'l = 2/3 - 100yi2/2, 

2/2 = ?/3 + 22/4- 1002/12/2 -2 -lO^yi, 

2/3 = -2/3 + 100yi2/2, 

2/4 = -2/4 + 10^2/i, 

t G [0,20], yi(0)= 2/2(0) = 1, y3(0) = 2/4(0) = 0, /iq = 2.5 • 10"^ 



The approximation of the Jacobian by a diagonal matrix is used when solving the test problems 
by AS0DE3. For the first test problem the diagonal matrix B with elements 611 = —0.013 — IOOO2/3, 
622 = —25002/3, 633 = — IOOO2/1 — 25OO2/2 are used. In the case of diagonal Jacobian approximation 
computational costs of additive methods are dominated by the number of right hand side function 
evaluations. So, computational costs of (2) per integration step are comparable to ones of explicit 
methods. Hence, AS0DE3 is compared with the following numerical codes based on well-known explicit 



RKM4 - 5-stage Merson method of order 4 [14], 
RKF5 - 6-stage Felberg method of order 5 [15], 
RKF7 - 13-stage Felberg method of order 7 [15], 

DPS - 13-stage Dormand and Prince method method of order 8 [16], 
and less well-known Runge-Kutta type method: 
RKN2 - 2-stage method of order 2 [12]. 

The overall computational costs (measured by the number of right hand side function evaluations 
over the integration interval) are given in the table 

Table 1. 

Computational costs of RKM4, RKF5, RKF7, DPS, RKN2, AS0DE3 with stability control. 





Tol 


RKA14 


RKF5 


RKF7 


DPS 


RKN2 


AS0DE3 


1 


10-2 


401 716 


401 005 


982 536 


717 526 


222 441 


9 351 


10-4 


400 627 


400 656 


982 150 


717 287 


222 481 


37 338 


2 


10-2 


13 391 594 


15 694 434 


38 429 196 


27 998 053 


8 682 849 


1 589 


10-4 


13 384 132 


15 691 105 


38 429 976 


27 993 793 


8 689 861 


7 711 


3 


10-2 


204 889 


237 942 


587 509 


431 591 


133 022 


3 129 


10-4 


206 647 


240 676 


565 396 


430 823 


132 987 


16 361 


4 


10-2 


10 832 


11 874 


29 991 


23 052 


6 585 


63 430 


10-^ 


10 236 


11 366 


28 819 


23 354 


7 627 


367 411 



8 Conclusions 

In addition to continuum mechanics problems, the constructed additive method can be used for solving 
locally unstable problems. In this case <f{y) corresponds to eigenvalues of the Jacobian matrix with 
positive real parts. In opposite to A-stable methods, explicit Runge Kutta methods are unstable in 
almost the entire right half plane and therefore are more suitable for detecting the local unstable 
solutions. For many locally unstable problems it is also easy to split the right hand side into stiff and 
non-stiff terms from physical considerations. 

So, in this paper, we constructed the third order additive method that is L-stable with respect 
to the implicit part and allows to use an arbitrary approximation of the Jacobian matrix without loss 
of accuracy. Automatic stepsize selection based on local error and stability control are performed and 
the auxiliary formulas for doing this were obtained without significant additional computational costs. 

The aim of numerical computations was to test the reliability and efficiency of the implemented 
integration algorithm with error and stability control as well as with diagonal Jacobian approximation. 
They didn't aim at solving practical problems of continuum mechanics and locally unstable problems. 
Numerical experiments show reliability and efhciency of the presented method. It follows from them 
that the method has good stability properties for solving mildly stiff problems and that the test 
problems turned out to be rather stiff for the explicit Runge-Kutta methods considered above. It is 
worth remarking that computational costs per step are comparable for both the additive method (with 
diagonal Jacobian approximation) and explicit ones. So, the implemented integration algorithm makes 
it possible to expend the range of applicability of explicit Runge-Kutta methods towards more stiff 
problems. 
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